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We study analytically the computational cost of the Generalised Hybrid Monte Carlo (GHMC) algorithm for 
free field theory. We calculate the autocorrelation functions of operators quadratic in the fields, and optimise 
the GHMC momentum mixing angle, the trajectory length, and the integration stepsize. We show that long 
trajectories are optimal for GHMC, and that standard HMC is much more efficient than algorithms based on 
the Second Order Langevin (L2MC) or Kramers Equation. We show that contrary to naive expectations HMC 
and h2MC have the same volume dependence, but their dynamical critical exponents are z = 1 and z — 3/2 
respectively. 



1. GENERALISED HMC 

The work reported here extends results pre- 
sented previously Jl|,||, to which the reader is 
referred for details. We begin by recalling that 
the generalised HMC algorithm is constructed 
from two kinds of update for a set of fields <f> and 
their conjugate momenta tt. 

Molecular Dynamics Monte Carlo: This 
consists of three parts: (1) MD: an approximate 
integration of Hamilton's equations on phase 
space which is exactly area-preserving and re- 
versible; U(t) : (<j), tt) h- > (<f>',Tr') where det U = 1 
and U(t) = U(-t)^ 1 . (2) A momentum flip 
F : tt i — > — tt. (3) MC: a Metropolis accept/reject 
test. 

Partial Momentum Refreshment: This 
mixes the Gaussian-distributed momenta n with 
Gaussian noise £: 
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The HMC algorithm || is the special case where 
i? = f . The L2MC/Kramers algorithm cor- 
responds to choosing arbitrary $ but MDMC tra- 
jectories of a single leapfrog integration step. 

1.1. Tunable Parameters 

The GHMC algorithm has three free param- 
eters, the trajectory length r, the momentum 
mixing angle d, and the integration step size St. 
These may be chosen arbitrarily without affect- 



ing the validity of the method, except for some 
special values for which the algorithm ceases to 
be ergodic. We may adjust these parameters to 
minimise the cost of a Monte Carlo computation, 
and the main goal of this work is to carry out this 
optimisation procedure for free field theory. 

Horowitz pointed out that the L2MC algorithm 
has the advantage of having a higher acceptance 
rate than HMC for a given step size, but he did 
not take in to account that it also requires a 
higher acceptance rate to get the same autocorre- 
lations because the trajectory is reversed at each 
MC rejection. It is not obvious a priori which of 
these effects dominates. 

The parameters r and i? may be chosen in- 
dependently from some distributions Pr(t) and 
Pm{&) for each trajectory. In the following we 
shall consider various choices for the momentum 
refreshment distribution Pm, but we shall always 
take a fixed value for §. 

2. FREE FIELD THEORY 

Consider a system of harmonic oscillators {4> p } 
for p £ Zy- The Hamiltonian on phase space is 
H = 5Ep e Z v (?l + ^W)- This describes free 
field theory in momentum space if the frequencies 
ujp are chosen as 
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3. AUTOCORRELATION FUNCTIONS 

3.1. Simple Markov Processes 

Let (4>i, 4>2, ■ ■ ■ , 4>n) be a sequence of field con- 
figurations generated by an equilibrated ergodic 
Markov process, and let (0(</>)) denote the expec- 
tation value of some connected operator O. We 
may define an unbiased estimator fl over the finite 
sequence of configurations by f2 = ■h Ylt=i ^(0t)> 

a i a ft n ti)\ (Wt+m^t)) 

As usual, we define Cn(<) = - — 7 \ — - as the 

(w) 

autocorrelation function for fi. The variance of 
the estimator is 
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where = Cn(£) is the integrated auto- 

correlation function for the operator f2 and A^ exp 
is the exponential autocorrelation time. This re- 
sult tells us that on average 1 + 2Aq correlated 
measurements are needed to reduce the variance 
by the same amount as a single truly independent 
measurement. 

3.2. Autocorrelation Functions for 
Quadratic Operators 

In order to carry out these calculations we make 
the simplifying assumption that the acceptance 
probability min(l, e~ SH ) for each trajectory may 
be replaced by its value averaged over phase space 
P acc = (min(l, e~ SH )); we neglect correlations 
between successive trajectories. Including such 
correlations leads to seemingly intractable com- 
plications. It is not obvious that our assumption 
corresponds to any systematic approximation ex- 
cept, of course, that it is valid when SH = 0. 

Details of the calculation of P acc for leapfrog 
integration are published elsewhere 



obtained by minimising the cost, that is by choos- 
ing f so as to satisfy dC/rff = dfL/dd = 0. 

We wish to compare the performance of the 
HMC, L2MC and GHMC algorithms for one di- 
mensional free field theory. To do this we com- 
pare the cost of generating a statistically indepen- 
dent measurement of the magnetic susceptibility 
M 2 , choosing the optimal values for the angle d 
and the average trajectory length f. We can min- 
imise the cost with respect to 1? without having 
to specify the form of the refresh distribution. 

The next step is to minimise the cost with re- 
spect to the average trajectory length £ = cut. 
Strictly speaking we should note that the accep- 
tance probability P acc is a function of f, but to 
a good approximation we may assume that P acc 
depends only upon the integration step size St 
except in the case of very short trajectories. 

4.1. Exponentially Distributed Trajectory 
Lengths 

To proceed further we need to choose a specific 
form for the momentum refresh distribution. In 
this section we will present results for the case 
of exponentially distributed trajectory lengths, 
Pr( t ) — re~ rT where the parameter r is just the 
inverse average trajectory length r — 1/f. 
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(-2P acc + l)c opt + (2P acc - l)c 
(^i-5P a „+4)c opt C Pt 2 
V + (-P .cc + 4)g opt 2 + (-4 + 3P acc )c opt 2 g opt 2 j 
P acc _<5™(P acc - l)c opt 3 £ opt 

' CoptP.co^TW - £opt_Pacc(5™C opt : 

-P acc <5™(-l + P acc )c opt £ opt 



4. COMPARISON OF COSTS 

If we make the reasonable assumption that the 
cost of the computation is proportional to the 
total fictitious (MD) time for which we have to 
integrate Hamilton's equations, then the cost (£. 
per independent configuration is proportional to 
(1 + 2Aq)t/St with f denoting the average length 
of a trajectory. The optimal trajectory length is 



This solution is a function of St and P acc which are 
not independent variables, and using the results 
for P acc (T, St) [0J|] we can compute the cost as a 
function of P acc as shown in Figure |l|. 

4.1.1. HMC 

The Hybrid Monte Carlo algorithm corre- 
sponds to setting d = tt/2, and we find that the 
optimal trajectory length in this case is £ opt = 
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Figure 1. Cost as a function of average Metropolis 
acceptance rate for the GHMC algorithm com- 
pared to HMC and L2MC for free field theory. 
The operator under consideration is the "mag- 
netic susceptibility" , i.e., the connected quadratic 
operator depending only on the lowest frequency 
mode. The corresponding parameters, the mo- 
mentum mixing angle $ op t and the average tra- 
jectory length measured as a fraction of the corre- 
lation length T opt /£ are also shown, all as a func- 
tion of the acceptance rate P acc . The inset graph 
shows the region where the acceptance rate is very 
close to unity which is where the L2MC algorithm 
has its minimum cost. 

1 / \A — -P acc , corresponding to a cost 



This is also shown in Figure ^. 

4.2. Fixed Length Trajectories 

For fixed length trajectories we shall only anal- 
yse the case of L2MC for which the trajectory 
length £ = ojSt. In this case the value of $ opt 
and the corresponding cost are also plotted in 
Figure [[} From this figure it is clear that the 
minimum cost occurs for P acc very close to unity, 



where the scaling variable x = V5t 6 is very small. 
We may then express c opt and -P acc as power se- 
ries in x, keeping only the first few terms. From 
these relations we find that the minimum cost for 
L2MC is 

Ct= {^y 4 v 5/4 m- 3/2 . 

This result tells us that not only does the tuned 
L2MC algorithm have a dynamical critical expo- 
nent z ~ 3/2, but also it has a volume depen- 
dence of exactly the same form as HMC 0,|. We 
may understand why this behaviour occurs rather 
than the naive l/ 7 / 6 m -1 by the following simple 
argument. 

If -P acc < 1 then the system will carry out a ran- 
dom walk backwards and forwards along a trajec- 
tory because the momentum, and thus the direc- 
tion of travel, must be reversed upon a Metropo- 
lis rejection. A simple minded analysis is that the 
average time between rejections must be 0(l/m) 
in order to achieve z = 1. This time is approxi- 
mately 

V P n (1 - P^nSr = ^t- = I. 

F or sm all 5t we have 1 — P acc = erf v 7 kV6r e cx 
VV5r e where k is a constant, and hence we must 
scale St so as to keep VSt 4 /m 2 fixed. Since the 
L2MC algorithm has a naive dynamical critical 
exponent z = 1, this means that the cost should 
vary as £ cx V^^mT 2 ) 1 ^ /mSr = V 5 / 4 m~ 3 / 2 . 
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